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Abstract. We consider mixing problems in the form of transient convection-diffusion equations 
with a velocity vector field with multiscale character and rough data. We assume that the velocity 
field has two scales, a coarse scale with slow spatial variation, which is responsible for advective 
transport and a fine scale with small amplitude that contributes to the mixing. For this problem 
we consider the estimation of filtered error quantities for solutions computed using a finite element 
method with symmetric stabilization. A posteriori error estimates and a priori error estimates 
are derived using the multiscale decomposition of the advective velocity to improve stability. All 
estimates are independent both of the Peclet number and of the regularity of the exact solution. 

1. Introduction. In spite of much progress in recent years the problem of deriv- 
ing a posteriori error estimates and a priori error estimates for transient convection- 
diffusion equations remain in its infancy. This is most likely due to the wide range of 
different problems covered by the equation class. Indeed depending on the character- 
istics of the velocity field and the molecular diffusion the solutions may feature very 
different behaviour. The difficulty will depend mainly on the smoothness of data and 
on the Peclet number 



where U is a characteristic velocity, L is a lengthscale and fi the molecular diffusion. If 
the variations of the transport velocity are small and the data are smooth, the solution 
will be smooth, with moderate Sobolev norm (at least the iJ ^norm) independent of 
the viscosity, and the flow will be computable using a standard Galerkin method 
for low Peclet number flows and stabilized finite element methods for high Peclet 
number flows. In this case L 2 -norm error estimates with an order can be established 
as we shall see below. Most difficult is the case of a high Peclet number and a 
strongly varying, or even turbulent, velocity field, transporting a concentration that 
is strongly fragmented and may dilute or concentrate. In this case a computation can 
experience strong amplification of errors due to repeated bifurcation of streamlines and 
an inexact representation of internal layers, due to spurious oscillations or numerical 
diffusion, leading to unphysical numerical mixing. This case is sometimes referred to 
as scalar turbulence and LES-models have been derived for the modelization of the 
passive scalar using filtering (see for example |18j). Such models encounter a similar 
Reynolds stress conundrum as the filtering of the full Navier-Stokes' equations, and 
hence the modeling error is very difficult to quantify. Another approach that has been 
attempted for this problem is heterogeneous multiscale methods (see for instance |14| ) 
however in that case the underlying theory is based on homogenization and completely 
depending on a periodicity hypothesis that in most applications will not hold. 

In this paper our approach is to apply a stabilized finite element method to the 
computation of the solution of the standard physical model, instead of a coarse grained 
model, the accuracy of the large scaled is measured by estimating the a regularized, 
or filtered error, equivalent to estimating the error in local averages of the solution. 
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The combination of these two ingredients allows us to derive error estimates with 
an order in h, for a norm that is in a certain sense in between H~ 1 and L 2 , but which 
contains the L 2 -norm of a filtered error. These estimates are robust in the sense that 
they do not depend on any high order Sobolev norm of the exact solution and they 
only have a exponential growth depending on the maximum gradient of the coarse 
scale velocity field, under a certain scale separation assumption given below. This 
means that the filtered quantities considered are robust under diverging fine scale 
characteristic trajectories. In some sense we extract the coarse scales for which we 
have some (provable) accuracy from the computation. 

The problem that we will consider takes the following form. Let fl be an open 
polygonal/polyhedral subset of R d , d = 2, 3 with boundary dfl, uq, / € L 2 (f2) and let 
e [C (0,T; W^°°(n))] d , such that V • = 0, ■ n dn \ dn = 0, fi > 0, then 

formally we may write, for t > find u <E Hq (fl) such that u(x, 0) = Uo(x) in ^ an d 

d t u + ■ Vu - fj,Au = /, in Q. (1.1) 

For the boundary conditions let u\gn = 0. We will denote the computational time 
interval by I := [0, T] and the space-time domain by Q := I x f2. The L 2 -scalar 
product over X, where X can be either a space domain of R d or a space-time domain, 
will be denoted by (-, -)x, the L 2 -scalar product over subsets X of will be denoted 
(•, -) x . In both cases the corresponding L 2 -norm is denoted by || • \\x- 

In this paper the analysis will be restricted to velocity fields with a particular mul- 
tiscale character. We will assume that € Vl /1 ' 00 (ri), such that V • = and that 
the problem is normalized so that U := ||/3||l~(q) = 1 and we assume that the char- 
acteristic lengthscale is given by L := 1, similarly we assume that T m L/\\0\\ l ^(q^ 
and hence also approimately one. Instead of making the standard assumption that 
||/3||w 1 .° o (0) is small, we assume that there is a decomposition of the velocity field, 

= + 0', 

where, for all t, \\0\\w 1 ' oo (n) ~ 1 and H/S'lHoom) ~ H where fi > is the molecular 
diffusion coefficient. This allows us to define a timescale for the flow relating to the 
coarse scale spatial variation and the fine scale amplitude, 

(tf)- 1 := \ supmmtpU-i), \\0'\\l~ { n)/ti ~ 1- (1-2) 

Essentially we assume that the velocity vectorfield can be decomposed in a coarse 
scale, responsible for transport, that is slowly varying in space and a fine scale, re- 
sponsible for mixing, that has small amplitude but may have very strong spatial 
variation. Expressed in Peclet numbers this means that the coarse scale Peclet num- 
ber may be arbitrarily high, whereas the fine scale Peclet number must be 0(1). A 
sharper value of tf given a molecular diffusion /i and a velocity field may be ob- 
tained by solving a certain minimzation problem in the L°°-norm that will be detailed 
in the a priori analysis. 

We will assume that the coarse scale velocity satisfies a pointwise non-penetration 
condition on the domain boundary, ■ uqq = 0. Here and in the following we used 
the notation a < b to denote a < Cb with C a moderate constant independent of the 
mesh-parameter and the physical parameters of the problem, (except those that are 
assumed to be unity.) We will also use the notation a ~ b for a < b and b < a. 

We will consider rough initial data or source terms, uq, f € i 2 (51), meaning that 
data can have strong fluctuation. This together with the high Peclet number and the 
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multiscale character of the velocity field may lead to complex, low regularity solutions. 
More precisely solutions are smooth, due to parabolic regularity, but with very large 
Sobolev norms, rendering standard a priori error estimates based on approximation 
theory worthless. Indeed classical global estimates for stabilized finite element meth- 
ods for time dependent convection-diffusion equations [HI [4J [8] yield the high Peclet 
number error estimate: 

|| ti - u h \\n + V(u - u h )\\ Q < Chi{l + Pe^ 1/2 )|w| L 2 (/;H 2 (0)) (1.3) 

where Uh denotes the approximate solution using piecewise affine approxomation. 
Even though |tt|flu(n) is huge in the presence of layers, stabilized methods are relevant 
for this case since one may derive localized error estimates, showing that perturbations 
can not spread too far upwind of crosswind in the stationary case [Ml HU [6] , or spread 
too far across characteristics in the transient case [IS]- The reason this works is that 
l M lff 2 (w) can be assumed to be small in a part of the domain w C O provided u has 
no layers in an h 5 -neighbourhood of u. In the estimates the bad part can be cut 
away using a suitably chosen weight function. This technique is not applicable in the 
present case, since the strong oscillations of the velocity field and the nonsmoothness 
of Mo and /, makes it unrealistic to assume that |it|#2( w is small in any part of the 
domain. The aim of this paper is to show that also in this case, stabilized finite element 
methods produce an improved solution compared to that of standard Galcrkin and 
that we can indeed derive error estimates for some large scale quantities defined by 
differential filtering. 

Drawing from earlier ideas on a posteriori error estimation (see [5J 1151 |3]) we 
propose to estimate a regularized error, or in other terms work in a norm in between 
H^ 1 and L 2 . Indeed let the regularized error e be defined by the partial differential 
equation 

-t)Ae + e = e(-,T), (1.4) 
with e\ga = 0, f) £ R + . We then prove that 

\^~4l + \\~e\\lf <C T {^- 



The constant Ct of the estimates is independent both of the Peclet number and of 
the regularity of the exact solution. The main contribution to Ct is a factor 

6XP (£) 

where T denotes the final time of the simulation. This factor limits the applicabil- 
ity of the analysis to time intervals of the size tf and hence limits the validity of 
the arguments to smooth vectorfields (3 or those satisfying the scale separation case 
discussed above. 

The parameter f) can be related to a filter width 5 by fj := S 2 , or an artificial 
diffusivity for elliptic smoothing. In both cases we see that the hidden constant 
includes dimensional quantities of order unity. For () constant the results can be 
interpreted as i? _1 -norm error estimates and in this norm the convergence rate is 
most likely optimal. 

Below we will first recall the standard L 2 -norm error analysis for high Peclet 
number flow and see what is required of data to obtain an estimate with an order 
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that is fully independent of /i (that is, we control suitable Sobolev norms of u a priori.) 
Then, we consider error estimation of filtered quantities and we obtain a posteriori 
error estimates. A priori estimates for rough solutions in the general case follow 
directly using standard stability estimates of the discrete solution. A similar analysis 
was carried out in the nonlinear case for the one dimensional Burgers' equation in 
[3]. In practice the a posteriori error estimates are likely the more useful, the upper 
bounds however show that the method always will converge. 

The weak formulation of equation writes, for t > 0, find u <E Hq{Vl) such 

that u(x, 0) = Uo(x) and 

(d t u,v) + a(u,v) = (/,«), WveH^(fl), (1.5) 

where (•, •) denotes the standard L 2 -product and a(-, ■) is defined by: 

a(u, v) := (j3 ■ Vu, v) + (/zVu, Vv). 

2. Stability and regularity of the solution. Since the constant in the es- 



timate (1.3) depends on the ff 2 -norm of the solution it is important to understand 
what quantities it is possible to control, both in general and in the special case that 
robustness in /i is required. In this section we will first show the standard global 
regularity estimate for parabolic problems detailing the dependence of the constant 
on /t. We will then show an estimate where the constant is independent of /i, but not 
of tf- In both cases we will need to assume that data have some smoothness. It is 
important to note that this regularity will not be required of the exact solution in the 
error analysis for filtered quantities that will follow. 



Lemma 2.1. (standard energy estimate) Letu be the solution of (1.1), then there 
holds 

Bup||u(t)|| n + ||m^Vu||q < / ||/(t)||ndfc + ||tio|| n . (2- 1 ) 
tei J i 



Proof. We take v = u in ( 1.5 ) and then use that 



(f,u) dt<sup|M,t)||n / \\f\\n dt 
i tei Ji 



from which the bound on sup tgJ ||u(£)||n follows. This bound is then used to get the 
if Abound, which is a consequence of the equation (1.5) and the equality 



||yLt 2 Vzi||g= / a(u,u) dt. 



□ 



Theorem 2.2. Let u be the weak solution of (1.1), with f G L (Q) and uq £ 



Hq (CI) assume that f2 is a convex domain, then there holds 



sup ||// 2 Vu(t)|| n + \\nu\\ L 2 {I . H 2 {n)) < [i 
tei 



< nr 1 / 2 



|/(t)||n#+||«o||n 



||/||Q + ||^ a Vu | 
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Proof. Multiplying (1.1) with — //Au, integrating over the space time domain 
Q* := (0, t*) x 57, with t* < T and applying the Cauchy-Schwarz inequality and the 
arithmetic-geometric inequality gives 



|At»V«(t*)lln + l|A*A«Ho. 



< ||^llioo { g. ) A*- 1 ||A**V«||^ 



Q* 



ImAu|| 



l/l 



Q* 



|/i3Vu ||n- 



It follows from (2.1| that 

|/3|Uoo (Q) /x- 1/2 ( / ||/(*)||ndf + ||«ol|n) + II/IIq + ||a*" V«o| 



sup ||/i2 Vii(t)|| n < 
te/ 



^AuHq < \\/3\\l~ { q)»- 1/2 ( / \\f(t)\\ nd t + \\u \\ n ) + ||/||q + \\^ V^lk 



and 



We conclude using elliptic regularity, recalling that the advection velocity is nor- 
malised. □ 

Remark 1. Observe that it follows that ||Vu||q < pT 1 ! 2 , sup tg j ||Vu||o < /i -1 and 
H u llL 2 (7;if 2 (f2)) /i -3 ^ 2 . Hence we conclude that the regularity estimates of Lemma 
\2.1\ and Theorem \2.2\ both are sensitive to the variation of the diffusivity and can only 
be used when Pc^ < 1. Also some regularity of the inital data is needed. Revisiting 
the estimate (1.3) in this context we get 

3/2 



\\pL*V(u-U h )\\ Q < 



( J ||/(t)|| n A+ Kiln + M 3 II/IIq + HmVuoIIo). 

(2-2) 

This estimate is clearly useless when Pe^ > 1. 

We will now show that using the scale separation of the advection field and assum- 
ing some more regularity of data, all inverse powers of the diffusivity may be avoided. 
The price to pay is an exponential constant, which however, under our assumption on 
the velocity field remains moderate. It should be noted that if no assumption is made 
on the velocity field the physical stability properties of the problem are very poor, 
since sharp, characteristic layers in the velocity may result in 0(/i5) width layers in 
the solution u and corresponding growth in the iJ -norm, making global estimates 
independent of \i impossible. 

Theorem 2.3. Let u be the weak solution of (1.1), with f G i?g(0) and uq G 
i?o(f2). Then for < f) there holds 

sup||nO,t)||„ + ||(b/i)*A U || Q +T-5||^Vn||Q + T-5||()^ tU ||Q 
tei 



with Ct — eV T -F/, where Tp is given by (1.2). 

Proof. Multiply equation ( 1.1 ) by — f)Au and integrate over Q 
t* < T, to obtain 



<C T (r 2 \\f\\Lm,m m) + I ||/||rft+||wo||f,), (2.3) 

[t*,T] xn with 



< 



(V/,J)Vu) c 



1 



f) 2 Vu | 



The second term in the left hand side does not have a sign and requires some further 
consideration. First split the velocity field in the large and the fine scale component, 



-(/3 • Vu, f)Au) Q . = - 09 ■ Vu, fjAu) g . - (/3' • Vu, f)Au) Q . , 

then integrate by parts in the term representing the large scale transport, noting that 
if t\ and i 2 denotes the two orthonormal tangential vectors to dfl, 

2 

(3 ■ Vu\ dn = (3 ■ n d n Vu ■ n d n\an + U(Vu ■ U)\ a n = 0. 

=o 1=1 =o 

-09 • Vu, f)A«)g. - (V(i9 • Vu), f,Vu) Q . 

Note that 

d d 

(V(0 - Vu), f)Vu) Q . - £((0 X4 j8) ■ Vu, f)d Xi u) Q . +209 - (d* 4 Vu), f^u) Q , . (2.4) 

i=l i=l 

For the first term of the right hand side we have 
d 

■ ^, W Xi u) Q , = ((V s j9) ■ Vu, t,Vu) Q , 

where Vs denotes the symmetric part of the gradient tensor. Similarly we obtain for 
the second part 

d d d 

^2(0- (d Xi Vu),t)d Xi u) Q * =^2^2(/3 j {d Xi d Xj u),t)d Xi u) Q , 

i—1 i—lj — 1 

d d d 

= ^2^2(Pj(d Xj d Xi u),l)d Xi u)Q^ = ^(/9- Vd Xi u,t)d Xi u) Q .. (2.5) 

i—1 

By the divergence theorem, recalling that (3 • n^a = 0, we have 

d ^ d 

^(/3 • Vd Xi u, \)d Xi u) Q * = --^(V ■ /3d Xi u, f)d Xi u) Q *. 
i=l i=l 

We conclude that, with X denoting the identity matrix, 

- 03 • Vu, f,Au) Q . = ((Vsi9 - l - V ■ jaZ)Vu, t)Vu) Q , (2.6) 

Observing that 

03' • Vu, f,Au) . < |/3'| M - 1/2 Vu|| Q , || (l)^Au\\ Q , 

<^l|f)'l/3'l«- 1/2 Vu||^+i||([)u)^Au||^ (2.7) 
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we have, 



r||^V«(-,t*)||^ + -||(Ai6)'A«||i.<& (V/.&Vu) 



where 



Defining 



ill^V^olln + (A(A/3',M)Vu, fjVu) (2.. 



A(/3,/3',/i) := V s /3- • (31 + *|/3'| 2 f\ 



T7 1 



1 



-sup _ inf a p (A(f3,(3',(i)) 

* tei /8£W I,DO (n) 



' /.• • — T -"4' »" '•:.< -mm. ,-' - }> • • (2-9) 

(3-n=o on an 

where <J p {-) denotes the largest positive eigenvalue of the matrix, we may write 



\\^Vu{;t*)\\h + \\m*M\Q* < t*\\l)^fV Q * 

+ ((f)- 1 + Tp^W^Vufg. + ||^V«o||^. (2-10) 

The claim regarding the control of the space derivatives now follows after an applica- 
tion of GronwalPs lemma, yielding 



\\t>^u(;t*)\\h + \\m^uf Q . < e 4 */^ (||^V/|||, + H^Vuoll 



(2.11) 



Combining this estimate with estimate (2.1) yields the claim for the first two terms 



in the left hand side of equation (2.3). Note that we also have 

2 < P T / f F 'n'.V n\2 



r-*||&*Vu||&<sup||&*Vu(-,t)||6 
tei 



(ll^V/|| 



t) 2 Vu | 



(2.12) 



For the control of the time derivative, simply multiply the equation with — tydtu and 
integrate to obtain 



l\\^d t uf Q < -{(3 ■ Vu, t)d t u) Q + h\tff\\% 



< \mL^(Q)\\^Vu\\ Q \\^d t u\\ Q + -\\0 f\\ 2 Q 



|f)*Vi*>| 



and we conclude using the estimate (2.12) and the observation that comparing the 
definitions (1.2) and (2.9) we have Tp < tf- □ 

Corollary 2.4. If in addition f2 is convex then there holds 



\{^) 2 u\ L 2 {I . H 2 {n)) < C T {h 2 \\f\\L 2 {I-H 1 (n)) 



l/ll*+ll«o||«i)- 



Proof. If is convex then elliptic regularity holds from which 

\{^)^u\ L 2 (I . H 2 [n)) < ||(/if))5Aii|| Q 

follows. □ 

Remark 2. It is straightforward to verify that under the assumption on (3 and (3' 
made in the introduction Tp ~ 1. 



3. Finite element discretizations. Let {Tk} be a family of nonoverlapping 
conforming, quasi uniform triangulations, Th '■= {K}h where the triangles K have 
diameter Hk and that is indexed by h := max/ijf. We let the set of interior faces 
{F}h of a triangulation Th be denoted by J ' . 

We will consider a standard finite element space of piecewise polynomial, contin- 
uous functions 

14° := {v h G : « h |jf G G %}, 

where Pk(K) denotes the polynomials of degree less than or equal to k on K. The 
following inverse inequalities are known to hold on Vh, 

\\Vv h \\ K < Cih^WvhWic (3.1) 

and 

\\vh\\dK < c t h~ 1/2 \\v h \\ K , . (3.2) 

The standard finite element method is then obtained by restricting the weak 



formulation (1.5| to the discrete space V®. For t > find Uh G V® such that Uh(x, 0) 
irhiio(x) and 

(d t u h ,v h ) +a(u h ,v h ) = (f,Vh), Vv h G V%. (3.3) 

Taking u/j = we immediately get the standard stability estimate 

i 

/ HAi'Vuhll^ dt) < f H/lln d*+||uo||n. 
Jo / Jo 

When the Peclet number is low this is a standard parabolic equation and it may 
be analysed using standard techniques, see Thomee [20] . 

However when the local Peclet number is high, the problem is a singularly per- 
turbed parabolic problem and the stability properties of the standard Galerkin method 
are in general insufficient for optimal convergence. In particular in the presence of 
layers the whole computational domain may be polluted by spurious oscillations, but 
as can be seen in Figure |3.1( convergence order is lost also for smooth solutions. 
In this case we study a Gaussian function convected one turn in a disc. Time dis- 
cretization is performed with the Crank-Nicolson method and we compare the result 
of the standard Galerkin method with those obtained using a symmetric stabilization 
method. 

In this paper we will consider symmetric stabilization methods only, however at 
least for time constant f3 one can obtain similar results for the SUPG-mcthod. 

3.1. Symmetric stabilization methods. In the last ten years there has been 
an important development in the field of high order symmetric stabilization methods. 
Such methods are typically obtained by the addition of a weakly consistent, dissipa- 
tive operator to the formulation. Some of the most important symmetric stabilization 
methods present in the literature are the subgrid viscosity method suggested by Guer- 
mond [H], the orthogonal subscale method proposed by Codina [9j, the local project 
method introduced by Becker and Braack [T] , the discontinuous Galerkin method [T7] 
and the continuous interior penalty (CIP) method suggested by Douglas and Dupont 
[10] and analysed by Burman and Hansbo [7]. The analysis that we propose herein 
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Fig. 3.1. Convergence in the L 2 -norm plotted against the timestep r = h, of the stabilized 
( dashed line for explicit and dotted line for implicit treatment of the stabilization operator) and un- 
stabilized (dash-dot) Crank- Nicolson method. The full line shows optimal second order convergence. 

will rely on an orthogonality argument and hence is valid only for the orthogonal 
subscales, the discontinuous Galcrkin method and for the CIP-method. To reduce 
the notational overhead we also exclude the DG-method from the discussion, but to 
partisans of this method the extension of the arguments should be straightforward. 
Herein we will assume that the following holds 

1. Sh(-,-) '■ Vh x Vh i-> R is a symmetric, bilinear and positive semi definite 
operator. 

2. For the i 2 -projection 7r^ : L 2 (il) i->- Vh there holds 

- Approximability 

\\h~ 1/2 (u - n h u)\\n + ||/iV(u - 7rfc«)||n < {h* + ^)h k \u\ H k+i {Q) (3.4) 

- Weak consistency of the stabilization operator 

s h (ir h u,ir h u) < c s h k+1/2 \u\ H k+i {n) (3.5) 

Enhanced continuity of the convection term: For all Vh <E Vh and /3 e 
W 1 '°°(Q), with ■ n dn \dn = 0, there holds 

\(u - TThU,f3 ■ Vu/,)| < ||/3|| w i,oo(n)||u - ir h u\\n\\vh\\n 

+ mi 00{Q) \\h- 1 / 2 (u-7r h u)\\as h (v h ,v h ). (3.6) 

As an example consider the CIP-method. Here s^(-, •) consists in a penalty on the 
jump of the gradient over element faces, and takes the form 

Sh(uh,Vh) ■= ^2 ( h F\\P- n F\\L™(F)Wuh-n F },lVv h -n F ]) F 

K FedK\Q 

where F denotes the faces in the mesh, [x] the jump of x over F, the orientation is 
not important, np a fixed but arbitrary normal associated to each face and (•, •) x the 
L 2 -scalar product over the d— 1-dimensional domain X. In the analysis below we will 
ufor simplicity use this stabilization operator. 



The stabilized finite element method then takes the form, for t > find Uh € Vh 
such that Uh(x, 0) = 717^0(0;) and 

(d t Uh,Vh) + a{uh,v h ) + s h (uh,v h ) = (f,Vh), Vt) ft . e 14. (3.7) 



For the satisfaction of the assumptions (3.4)-(3.6) we refer to [5]. Although in that 



reference Nitsche-type boundary conditions are used, the same argument works when- 
ever (3 ■ V^/Jan — 0, which will the case herein. For completeness we show how to 



modify the result of [5j to obtain (3.6| in our case. We also show that the flow field 
(3 in the stabilization operator may be replaced by (3 at the cost of a nonessential 
perturbation. The key result is the following Lemma. 

Lemma 3.1. Assume that (3 ■ nan | an = then there holds 

inf ||/ii(7r /3 • Vu h - Vh)\\n < s h (u h ,u h )i + hi\\j3\\ W i,^(n)\\uh\\n, 
v h ev h 

where ttq/3 is a piecewise constant approximation of (3 that will be defined below. 

Proof. Let wof3 be the projection onto element wise constants such that for every 
K that has no face on the boundary there holds 



7To/3 dx = I (3 dx. 
k J K 

On elements adjacent to the boundary define ttq(3 by 



tto/3 • non dx = (3 ■ ng n dx 

dKudQ JdKUdn 



and 



/ 7To/3 • U dx = / /3 ■ tj dx, i = 1, d — 1. 

JdKudn JdKudn 

It follows from standard approximation that for all K £ Th, 

11/3 - koP\\l°°(K) < WPWwi^iKjh-K- 
Note that for any element K with one face on the boundary there holds 

2 

Tr f3 ■ Vuh\au = tto/3 • ndnVu h ■ n d n\aa + ^7r /3 • tj( Vuh ■ U) \dn = 0. 
=0 1=1 =0 

It then follows using the same arguments as in jS] that 

inf \\hi{TT l3-Vu h -v h )\\n< I V \\hliT p ■ Vu h jf F ) . 

By adding and subtracting (3 inside the jump and by applying a triangle inequality 
followed by a trace inequality we have 

I WhfroP- Vu h ]\\ 2 F I < s h (u h ,u h )2 +hi\\P\\ w i,~ i n ) \\u h \\a 

\F£T ) 
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and the proof is finished. □ 

The continuity (3.6) is now obtained by adding and subtracting n f3 in the right slot 
of the left hand side of the equation 

\(u-ir h u,/3-Vv h )\ < \(u-Tr h u,(/3-TT /3)-Wv h )\+ inf \(u-ir h u,ir [3-Vv h -w h )\ 

WhEVh 

< ||^|| w i,o«.( n ) ||u - Tr h u\\n\\v h \\ n 

+ \\P\\L^(n)\\h~ 1/2 (u - Tfhu)\\nSh(vh,Vh)- (3.8) 

where a Cauchy-Schwarz inequality, an inverse inequality, the approximation proper- 
ties of 7To/3 and Lemma |3 . 1 1 have been used in the last inequality. 
Lemma 3.2. Let 

Sh(uh,Vh) ■= X! ( h F\\P • nF\\L°°(F)$yuh ■ n F ], [Vv h ■ n F j) F . 

K FeJF 

Then 

s h (u hl u h ) < s h (u h ,u h ) + h^\\fi^Vuh\\ci 



s h (u h ,u h ) < s h (u h ,u h ) + h 2 \\n 2 Vu h \ 



2 



Proof. The proof for the upper and the lower bounds are similar, so we consider 
only the first inequality. By using the decomposition of (3 and a triangular inequality 
we have 

Sh{u h ,u h ) < 7 ^2 II Ml/ 3 ' MU~(F)[Vu/i]||| +7 ^2 1 1 M 1/3' ' n F \\ L oc^a ) [Vu h ]\\%. 

Using now the size constraint on and a trace inequality we conclude 

7^ || Ml/3'- n F \\ L aa ln) \yu h ] III < 7 W h K^u h f K <-fhi\\^Vu h \\ 2 Q . 

F6J KeT h 

□ 

3.1.1. Stability and convergence results. We will now recall the main results 
on stability and convergence of symmetric stabilization methods. These results are 
minor modifications of those in 0]. For convenience we introduce a triple norm 
associated to the stabilized method. Let 



u h\\\h : = / (iM Vu^ll 2 + s h (u h ,u h )J dt. 



Following the proof of Lemma [2T| it is straightforward to derive the following stability 
and error estimates. 

Lemma 3.3. (Stability) 
Let Uh be the solution of (3.7), with 7 > and jbc > Co > 0, cq large enough, then 
there holds 

swp\\u h {t)\\ n + \\\u h \\\ h < / ||/||n dt + |K||n- 
tei Jo 
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Applying stability to the error equation leads to the following error estimate: 

Theorem 3.4. (Convergence CIP-method) 
Let u G L 2 (I;H r (n)) be the solution of (1.5) and Uh the solution of (3.7). Then 
there holds 

sup||u-Ufc(t)|| n + |||u-Ufo||U < {Tp X T?h+(Tp*h? + l)fti + Ai^^M^/^tfi)), 
tei 

with s — min(r, k+ 1) and where we used that ||j9||.l°°(q) ~ 1- 

Proof. For simplicity we only give the proof in the form of a final time result. Let 
c-h '■= Uh — 7T/jM, using the same arguments as for the stability Lemma |2.1| we have 

\\\ e h{T)\\ 2 n + \\\e h \\\l = J[(d t e h ,e h ) n + a(e h ,e h ) + s h (e h ,e h )} dt. 

Using Galerkin orthogonality and the orthogonality of the i 2 -projection we have 



l\\e h (T)\\l + \\\e h \\\l 



~a(u - n h u, e h ) + s h (-K h u, e h )] dt. 



Using the decomposition of the velocity we have 

a(u-n h u,e h ) < (u- n h u, (3 ■ Ve h ) + (u- -K h u, ■ Ve h ) + \\^ V(n- 7r ft u)||n||At' Ve^|| fi 



and using the continuity (3.6) in the first term of the right hand side and the bound 
on in the second we have 

a(u- ir h u,e h ) < ||V/3||ioo(n)||(u - 7r h u)||n||e/i|| n 

i 1 



nhu\\n + Ha* 2 V(u - 7rhu)||n)|||e/i|||h. 



It follows after a Cauchy-Schwarz inequality and an arithmetic-geometric inequality 
in the right hand side that 

\\e h (T)f Q + \\\e h \\\l < \\VP\\l™ {Q) T\\h- 1/2 (u - n h u)\\ 2 Q 

+ (Tp 1 + WPh-^h-^Wu - n h u\\ 2 Q + \\^V{u - n h u)\\ 2 Q + T- l \\e h \\ 2 Q . 



We conclude by applying Gronwall's lemma and the approximation results of (3.4). □ 
For high mesh Peclet numbers this estimate is optimal in the L°°(I; L 2 (f2))-norm and 
for low mesh Peclet numbers it is optimal in the L 2 (I; 7J 1 (r2)-norm. In the latter 
case the convergence in the L 2 -norm can be improved under certain assumptions on 
the time variation of f3, see [5D]. Note that this estimate does not have exponential 
growth, however that factor is hidden in the Sobolev norm of the exact solution as we 



shall see. Combining this convergence result with the regularity result of Theorem 2.3 



we may prove the following estimate that is fully independent of /i, in the sense that 
we also control the Sobolev norm in the constant. Note however that smoothness of 
the source term and the initial data is required. 

Corollary 3.5. Let f e L 2 (I; Hq(Q,)) andu € Hq(Q), letu denote the solution 
of (1.5) and Uh the solution of (3.7) and assume that Pe^ >> 1. Then 



\\(u-u h )(T)\\ n < CriTHl + Tp 1 ) + l)hH\\f\\ L 2 (I . mm + \\Vu \\n). 



Proof. An immediate consequence of the Theorem 3.4 and Theorem 2.3 □ 
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4. Perturbation equation and the dual problem. The analysis uses a per- 
turbation equation and an associated dual problem. In this work we use the pertur- 
bation equation associated to the weak form of the equations (1.5| and the standard 
Galerkin formulation (3.3 1, for the perturbation equation. Taking the difference of 
the two formulations, setting e — u — Uh and integrating by parts we obtain 

(d t e, if) + o(e, if) = (e(T),ip(T)) - (e(0), <p{Q) - (e, d t <p + /3- Vtp) + ( M Ve, V<p). 

This suggests the adjoint equation, find ip € Hq(Q) such that 

—dtf — (3 ■ — fiAip — in f2 

ip = Oonffl (4.1) 
<p(-,T) = ttQinfi. 

So that the following error representation holds 

(e(T), *) = (e(0), ip(0) + (d t e, ip) + a(e, ip). (4.2) 
Observe that since e\gn = and <p\dn = we have dt^p\on = 0, 

(f3 ■ Vip)\en = P ■ ngn Vtp • n dn +/3x n m ■ yip x ngn = 0. 

=o =o 

Finally, by the equation we have A^l^n = 0. These properties will be useful to prove 
the necessary bounds on the dual solution. The properties of the continuous dual 
problem will be crucial for our argument. We will now proceed and discuss the choice 
of <!' and the associated stability estimates on ip. 

4.1. Regularization of the error and weak norms. Since it appears not to 
be possible to prove a posteriori error estimates for the L 2 -norm independently of the 
Peclet number, unless one ressorts to a saturation assumption, we will here consider 
a regularized error, where a parameter t) that may ultimately depend on h sets the 
scale of the regularization. We recall the problem (1.4 1 for a given computational 
error e := u — Uh, find e such that e\ga = and 

-f)Ae + e = e(-,T). 

On weak form the problem writes: find e € Hq(£1) such that 

(f)Ve,V«) + (e»=(e,i>), V«e£#(fi). (4.3) 

This is what is commonly called a Helmholtz filter or a differential filter, although it 
is not properly speaking a filter. The key observation here is that when t) is small, 
the filtered solution is close to the solution where ever the solution is smooth. Close 
to layers or other strongly localised features of e, u — u may be large locally, also for 
f) small. Associated with the problem (4.3) we have the norm 



\~e\\l 



= l|f) 5 Ve|| 2 



lell 2 



and an associated relation, obtained by testing (4.3) with v = e, 

(e,e). 



(4.4) 



We deduce from (4.4) that choosing ^ 
representation for ||e||?: 



e in (4.1 ) above leads to the following error 

(4.5) 



(e(0),<p{0)) + {d t e,ip)+a(e,<p). 
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4.2. Stability of the dual solution. The advantage of using the dual technique 
is that instead of relying on regularity estimates for the exact solutions we can use 
regularity of the adjoint equation, which may be better behaved provided the data is 
well chosen. The following Theorem gives a precise characterization of the regularity 
of the dual problem in the multiscale framework. 



Theorem 4.1. Let tp(x,t) be the weak solution of (4.2), with 4* = e then there 
holds 



sup||^V^(-,t)||n +T- 1 ||f)3V^|| Q + T- 1 \\t ) 1 id t <p\\ Q + IK^OWllo £ Cr||e|| b , 
tei 



with Ct = e^ Tp ' where Tp is given by (1.2). If in addition ft is convex, then 



\(i)^)^Lp\ L 2 {I . H 2 (n)) < C T \\e\\ 



Proof. The dual problem is equivalent to the forward problem after a change of 
variable t = —t and x = —x, with the source term / = and the initial data uq = e. 
The result then follows from Theorem |2.3| and Corollary |2.4| □ 



5. Error estimates. In this section we will prove estimates for \\e\\t, where the 
constant is robust in fi (under our assumptions on the data). We will only consider 
the case of semi discretization in space and show how to prove a posteriori error 
estimates, where the stability constant is essentially Ct of Theorem [44] The a priori 
error estimates then follows using the fact that the a posteriori residuals are a priori 
controlled by the discrete stability estimate (3.3). We will first consider the case 



where no convexity assumption is made on the domain so that elliptic regularity can 
not be assumed. Then we will consider the case of insufficient data, i.e. when only 
/3 is known, and show that under our assumptions on data we can obtain an upper 
bound of the error also in this case, where a nonconsistent part limits the asymptotic 
convergence, in the regime that we are interested in however this part is smaller than 
the discretization error. 

In the low Peclet number regime, we need to modify our estimate to obtain 
optimality. These modifications are possible only if we assume that the domain is 
convex so that elliptic regularity holds. In that case the estimates below can be made 
robust also in the limit of low Peclet numbers by straightforward modifications. We 
omit the details. 



5.1. A posteriori and a priori error estimates. To prove a posteriori error 
estimates we use a duality technique together with the a priori control of the dual 
solution. In practice one may need to ressort to numerical solution of the dual problem. 
In the multiscale framework it is however unclear if the norm of the gradient of the 
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dual solution is computable or not in general. 

Theorem 5.1. (A posteriori error estimate) Let e be defined by (4.3 1. Then 
there holds 

\\eh<C T (^) 2 ( ! inf \\h2(p-Vu h -v h )\\ n dt 

+ I I inf Yl WhHftAuh - v h )f K ) dt+ f ( V WnlVuh-nplfp ) dt 
J l\ VheVh KeT h J Jl \Fer J 

+ J s h {u h ,u h )^dt + hi J \\f-TT h f\\u dt + h?\\uo-n h u \\ n y (5.1) 



Proof. Starting from (4.5 ) and using Galerkin orthogonality and the orthogonality 
of the L 2 -projection we have 

= (e(0), ip(0) - n h (p)n + (d t e, ip - ir h <p) Q + J (a(e, <p - n h <p) + s h {u h , n h <p)) dt. 

Using the weak formulation (1.5) and the orthogonality of the L 2 -projection we may 
write 

= (/ - Khf, KhW - <p)q + (e(0), <p(0) - Kh<p(0))n + (/3 ■ Vuj, - v h ,w h (p - tp) Q 

+ (pVuh, V(7r/ l y- <p))q + J s h (u h ,ir h tp) dt. (5.3) 

Considering the right hand side term by term we get using Cauchy-Schwarz inequality, 
approximation and the stability of Theorem |4.1| 



(f-rt h f,ir h tp-tp) Q <hf) 1/2 / ||/ - TTh/Hndi sup ||fp Vy>||j 

J i tei 



~. ■('■; ( ~Y ^Wf-ithfWndt ||e||„, (5.4) 



(e(0), <^(0) - ir h <p(0))a < hi)- 1/2 ||e(0)|| n sup ||^V^||n 

tei 



{/3-Vu h -v h ,Tr h <p-<p)Q <h*\) 1/2 / \\h2([3-Vu h -v h )\\ n disup||fpV¥>(',*)||n 

J i tei 
i 

S^lrT f \\hHP-Vu h -v h )\\ u dtp||„. (5.6) 



i 



In the fourth term we first integrate by parts on each element and then proceed with 
trace inequalities, followed by approximation for the dual solution, 



{nVu hl V{-K h if - <p)) Q 

< [ ( mf ]T \\hH»Au h -v h )\\ 2 K ) + I \HVu h -n F }f F ] dt 



sup hhf - vIIf + \\h- 1/2 {ir h <p - ip)\\f, dt 



SGI 



\ F £ F 

i 

I R(u h ,fi)Bup\\tfV<p(-,t)\\n 



<Ct^J n{u h ,ix)\\eh. (5.7) 

Finally the stabilization term is handled using a Cauchy-Schwarz inequality, followed 
by a trace inequality and the iZ^-stability of the L 2 -projection on quasi-uniform 
meshes. 



' s h (u h ,Tr h ip) dt < \\P\\lao (0) / s(u h ,u h ) 2 dt h 2 sup || V(/3(- , t) || a 
/ ' J i tei 



: C T I j-J J^u^Uh)* dt\\~e\\ h . (5.8) 



The claim follows by collecting the upper bounds ( 5.4 1-( 5.8 ) and dividing by |je||p). □ 
Theorem 5.2. (A priori error estimate) Assume that Pe/j > 1 then there holds 



MI*;$Or(rJ 0- + h * + T5 ) [J ll/lln*+Wn 



Proof. The result follows from estimate (5.1) by bounding all the residual terms 



using (3.3). First observe that since the Peclet number is high we may use inverse 



inequalities and trace inequalities to show that 



n{u h ,n)<T*\\^u h \\ Q <T* I J ||/||ndi + ||«o||n 



For the contributions on the faces we have used 



jf \HVu h ■ n F ]\\ r dt < Ml^h^ ^ £ {Vu h ■ n F j\\ 2 F dt 



< 



||/3|IS- (n) ||/i»V« fc ||o<|||ttfc||| h . (5.9) 
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Similarly using a Cauchy-Schwarz inequality in time we have 



s(u h ,u h )i dt<Ti I s(u h ,u h ) dt<Ti I / ||/|| n dt + ||uo||n 



We finally consider the first term on the right hand side of (5.1 1. First note that 



/ inf ||^(/3.Vti h -t; fc )||ndt<Ti inf \\h^{f3-Vu h -v h )\\ Q +T^\\h^ -Vu h \\ Q . 
J i v h ev% v h ev h 



By Lemma 3.1 and Lemma 3.2 we have 



inf \\hi((3-Vu h -Vh)\\Q<h 2 \\l3\\ W i,eof n )\\uh\\Q+[ s h (u h ,u h )dt 
v h ev h ■ \jj j 

< h'\\l3\\ w i,oo(Q)\\uh\\Q + h^\\^Vu h \\ Q + (J s h (u h ,u h )dtj 

It follows that 

inf \\hi(0-Vu h -v h )\\ Q 

Vh£V h 

< mtn(hiTl\\0\\ w i.ao [a) , ||/3||^ (Q) )(sup |K(-,t)|| fi + IKHU)- (5.10) 
Using the assumption on the small scale fluctuations ||/3'|||oo(q) < f- we have 



We conclude by collecting terms and applying Lemma |3.3| □ 

Remark 3. (The necessity of stabilization for robustness) Note that the stability 
of the dual problem holds regardless of the numerical method used. The stabilization 
in the numerical method allows us to control the first residual in the a posteriori 



error estimate, by using the discrete stability estimate (3.3). If no stabilization is 
present there is no control of the streamline derivative making it impossible to obtain 
uniformity in fi. Another observation that is worthwhile is that the above argument 
is valid only for high macroscopic Peclet number. One may prove that the situation 
improves if the domain is convex so that elliptic regularity can be used, in particular 
we then get an estimate that is valid also in the low Reynolds number regime, but also 
in this case stabilization is required to yield uniformity in fi. 

5.1.1. The degenerate case of unknown . In many relevant cases may 
be unknown, or only statistics of it known. Then some stochastic method may be 
used to recover expectancy values for the solution. In this section we will consider 
the situation, that is simply excluded from the computation and we will show that 
under our assumptions on the small scale velocity fluctuations the approximation 
result still holds. Indeed in the high Peclet number regime the consistency error 
made by dropping the fine scale fluctuations of (3 is smaller than the discretization 
error. 

Here we use an advective field (3 that we assume is divergence free and let (3 be 



replaced by (3 in (3.7) 
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Theorem 5.3. Let u be the solution of (1.5) and Uh the solution of (3.7), with 
(3 instead of [3 for the advective field, assume that Pe/, >> 1, then 



\\e\U <C T 



\\f\\u dt + \\u \\n 



) 



Proof. We proceed as in the proofs of Theorems |5.1| and |5.2| 



e 



2 _ 



(e(0), (p(0) - TT h ip) n + (d t e, <p - ■kh&q 



i 



;). (5.11) 



The only thing that differs from the previous analysis is the last term in the right 
hand side. Except for this term the proof proceeds as before. We will therefore here 
only show how this term may be bounded. After an integration by parts we have 



(0 ■ Vu h , tp) Q = (u h ,f3' ■ V<p) Q < sup \\u h (;t)\\ n T\\f3'\\ L ^ {Q) i)- 1/2 sup \\tfv<p(;t)\\n 



It follows that the consistency error is of the same order as the discretization error. 
If the Peclet number is large, the contribution from the discretization error can be 
assumed to be dominating and the same order of convergence as in the unperturbed 
case should be observed, until the Peclet number becomes so small that the inconsis- 
tency dominates. This means that in the high Peclet regime, if data are known to be 
rough, noise in the velocities satisying the constraint on may be neglected. Hence 
in any stochastic analysis or computation at high Reynolds number fluctutations must 
be larger than /Z5 to be essential. 

5.2. Conclusion. We have derived robust a posteriori and a priori error esti- 
mates for transient convection-diffusion equations. The upshot is that the estimates 
are completely robust with respect to the Peclet number, in the sense that we also 
control the Sobolev norms of the exact solution in the error constant. The estimates 
allow for low regularity data and multiscale advection, that may have strong spatial 
variation on the fine scale under a special scale separation condition. It is easy to see 
that a general high Peclet transport problem with very strong spatial variation in the 
advection velocity must be ill-posed since a small perturbation of the initial position 
of a flow particle may have an important influence on its trajectory. So our work is 
an attempt to separate computable cases from inherently unstable ones, using scale 
separation. In this spirit, our paper proposes a first step towards an understanding 
of what transport problems are computable in the high Peclet regime, beyond the 
standard assumption of smooth data and sets a baseline for what should be achieved 
theoretically in the analysis of more involved methods, such as multiscale methods, 
in order to claim that they produce an accuracy beyond what is obtained using a 
standard stabilized finite element method. 



Note that under our assumptions ||/3'||l°°(q) 



i - i 

< \\0\i^(Q) h * leading to 




□ 
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